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Abstract 

We demonstrate how the mixed dynamic form factor (MDFF) can be interpreted as a quadratic form. This makes it 
possible to use matrix diagonalization methods to reduce the number of terms that need to be taken into account when 
calculating the inelastic scattering of electrons in a crystal. It also leads in a natural way to a new basis that helps 
elucidate the underlying physics. The new method is applied to several cases to show its versatility. In particular, 
predictions are made for directly imaging atomic orbitals in crystals. 
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1. Introduction 

Nowadays, simulations are indispensable both for plan- 
ning and for interpreting experiments in the transmission 
electron microscope (TEM), in particular when working 
with electron energy loss spectrometry (EELS). The key 
quantity for simulating inelastic electron scattering is the 
mixed dynamic form factor (MDFF) [UJ3]. In many cases, 
this complex quantity is simplified by several approxima- 
tions, like, for instance, the dipole approximation. Re- 
cently, it has been shown, however, that this can lead to 
quite severe errors [4]. Furthermore, with recent advances 
of aberration corrected microscopes, more accurate cal- 
culations of the MDFF will become essential for future 
experiments. 

In this work, we will give a brief repetition of the mixed 
dynamic form factor, followed by an analysis of how a ba- 
sis transformation can bring it into a simpler form that 
is much better to handle numerically. Furthermore, the 
physical significance of this procedure will be outlined. In 
the last part, the new formalism will be applied to both 
existing and new measurement setups to study its appli- 
cability and versatility. 



2. The mixed dynamic form factor and its factor- 
ization 

In the most general approach, the quantum mechanical 
system consisting of both the probe electron and the sam- 
ple can best be described by a density operator p or its ma- 
trix elements, the so-called density matrix p [5]. Adopting 
the density matrix formalism instead of the simpler wave 



function approach is greatly beneficial as one cannot ob- 
serve the target's final state directly. This ignorance of a 
part of the system after an inelastic interaction gives rise 
to a mixed state which can be described very effectively 
using the density matrix [2j O [5] . 

Before the interaction, the probe and the target system 
can be considered independent. For the sake of simplicity, 
we will furthermore assume that both systems are initially 
in a pure state, i.e., each can be described by a single wave 
function. Then, the density operator of the whole system 
before the interaction 
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where (g) denotes the direct product. Throughout this arti- 
cle, we use small letters when referring to the probe beam 
and capital letters when referring to the target. 

In first order Born approximation, the density operator 
after the inelastic interaction mediated by an interaction 
potential V is given by 
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Since the target system is not observed directly, one has to 
construct the reduced density operator for the probe beam 
by summing incoherently over all possible final states of 
the target. This reduced density operator is given by 



p = ^2(F\V\I)\i) (i\(I\V^\F) 



(3) 
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which can then be propagated elastically through the crys- 
tal and be used to predict the outcome of measurements 
in different geometries. It must be emphasized that the 
ordering of the terms is vital here, since V in general acts 
on both the probe and the target states, which results in 
an entanglement of the two. 
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In EELS experiments, the interaction operator V is the 
Coulomb interaction operator. Its two most common basis 
representations are in configuration space, 

V(r, r') = (r\V\r f ) = S ^~^h {r - r 7 )*^/ - E F + E) 
\r-R\ 

=: V(r)5(r-r'), 



and in reciprocal space, 
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bound initial states that give rise to EELS core losses. 
Hence, the initial state written as |/|jj^)R while the final 
states are expanded in terms of \LM^S). In the follow- 
ing, we will also sum incoherently over j z since that quan- 
tum number of the initial state is typically unknown. In 
the Kohn-Sham approximation, the MDFF is then given 

by [SHE] 

SW)=£ E (F\LM l -S){LM l -S\V{q)\l\ 

Fj z LMSL'M'S' 
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The constant prefactors were omitted for the sake of clar- 
ity. Ei,E F denote the energies of the target's initial and 
final states, respectively, and E is the energy loss of the 
probe electron. 

In these two representations, the reduced density op- 
erator reads 

p(r,r')=E / dr(F\(r\V\f)\I) (f\i) 

F J 

Jdf (i\r') (I\ (f'\V^\r')\F) 
= J2(F\V(r)\I) (I\VHr')\F) (r\i) (i\r') 

F 

=S(r,r') (r\i) (i\r') 
p(k, k') = E / dfe (F\ (k\V\k) \I) (k\i) ( 6 ) 
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J dk (i\k) (I\(k\V^\k f )\F) 

=E JJ dqdq ' w^ww 

(k + q\i) (i\k' + q') 
= U dqdq f S(q, q f ) (k + q\i) (i\k' + q') . 

Here, the MDFF S(q, q') and the real-space MDFF (rMDFF) 
S(r,r f ) were introduced which are related by a Fourier 
transformation^] It is noteworthy that — due to the par- 
ticular properties of the Coulomb operator — the rMDFF 
can be multiplied on the initial probe wave functions, whereas 
the MDFF has to be convolved with them. 

In order to perform calculations, one not only has to 
specify a basis for the probe states, but also for the tar- 
get states. Usually, one chooses a spherical harmonics ba- 
sis which is particularly useful for describing the tightly 
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Here, 



(j\(q))ELSj = J uf s (r)jx(qr)Rij(r)r 2 dr 
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is the weighted radial wave function overlap [6j [9] with the 
initial state's radial wave function Rj(r), the final state's 
radial wave function uf s (r) and the spherical Bessel func- 
tion j x . The Y.kn D lms ( d l?M'S'Y is the cross density of 
states (XDOS) and the ^* * *^ are Wigner 3j symbols. 

While this choice of basis is very convenient as a start- 
ing point (as it is used, e.g., in WIEN2k [E]), it is by no 
means the only or the optimal choice. This can be seen by 
collecting terms depending on q and terms depending on 



2 This takes into account the spin-orbit coupling of the tightly 
bound core states 



1 Contrary to the convention adopted in previous works, we in- 
clude the l/q 2 q' 2 term in the definition of the MDFF as it makes 
the definition more concise and easy to use. 
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q' . With the abbreviations 

a =(A, /i, L, 5) 
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the MDFF can be rewritten as 



olol' 

= 9 (q)t.S.g(q'), 



(9) 
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where the matrix S collects all qr, q' independent terms 
and can be computed in a straight-forward way once the 
XDOS is known (e.g., from DFT calculations). The g in 
turn can be interpreted as a vector of functions. This is a 
well-known quadratic form. In particular, S is hermitian 



(as is shown in Appendix A). 

With the default settings, WIEN2k produces data with 
< L, V < 3. When including transitions up to quadrupole 
order (A = 2), S is a 72 x 72 matrix, resulting in up to 
5184 terms in the MDFF that in principle would all have 
to be handled separately. In practice, some of the entries 
vanish due to selection rules, while for some others the 
hermiticity of S can be exploited. Still, many off-diagonal 
elements generally remain. These off-diagonal elements 
imply correlations between the basis vectors [3] and hence 
represent additional information (e.g., symmetries) about 
the underlying system that can be used to simplify the 
problem. 

To exploit this additional information, one can insert 
unitary matrices U in the following way: 

S(q, q') = fl(g) f • U^U ■ S • U*U ■ g(q'). (11) 

Since for any hermitian matrix, a unitarian matrix exists 
such that USU^ is a diagonal matrix .D, one only has to 
find such a U. This is straight forward using, e.g., eigen- 
value solvers, a singular value decomposition, or a Schur 
decomposition. With the abbreviation g(q) = U-g(q), 
the MDFF becomes 



In terms of quadratic forms, the transformation U is 
a principal axis transformation. In quantum mechanical 
terms, it is a basis transformation into the eigensystem of 
the MDFF. In essence, it recovers the "physical" basis of 
independent — i.e., uncorrelated because of vanishing off- 
diagonal terms — transitions^] With the default settings 
of WIEN2k, this means that the problem was reduced from 
at most 5184 terms to at most 72 terms. 

In numerical simulations, it is usually beneficial to work 
with the rMDFF as it can be multiplied directly onto the 
incident density matrix. Since the rMDFF is related to 
the MDFF by a Fourier transformation and the q and 
q' dependencies have been decoupled, the rMDFF simply 
reads 

S(r,r , ) = ~g(r)^D-~g(r') (13) 

with the same matrix D as for the MDFF and g(r) = 
FT q [g(q)]. By renormalizing the g such that g a (r) := 
VD aa g a (r), the MDFF can be further simplified to 

S(r, r') = » (r)t • g(r') = £ fl a (r)*fl a (r') (14) 

a 

Hence, the reduced density matrix of the probe electron 
after the inelastic interaction in configuration space can be 
written as 



(15) 



where we wrote <p(r) = (r\i) for the incident probe electron 
wave function. It is quite obvious that the diagonalization 
of the rMDFF has resulted in a factorization of the density 
operator p = ^ a \a) (a\. Finally — when measuring a real 
space image — , the measurable intensity I is given 



I(r) = p(r,r) = Y / \9a(r)<p(r))\' 



(16) 



In the description above, elastic scattering of the probe 
beam after the inelastic scattering event has not been in- 
cluded for the sake of simplicity. Since each \a) in itself is 
a pure state, it can be propagated elastically through the 
rest of the crystal with existing methods (e.g., the mul- 
tislice approach [H]), however. While this does change 
the image, of course, it does not change the underlying 
physics, or the coherence properties and the general state- 
ments made above still hold. 



S(q,q')=g(qy-D.g(q'). 



(12) 



3 Note that, depending on the EELS-edge and multipole orders 
considered, it may or may not be possible to determine the tar- 
get's "physical" basis \F) from these transitions. Considering, e.g., a 
dipole-allowed transition from an initial p state to a final d state, one 
has only three transition elements (/x G { — 1, 0, 1}), but 5 final states. 
Under these circumstances, not all information about the final states 
can be probed, unless one takes into account other multipole orders 
or final states. 

4 Here, an ideal lens system was assumed. Real lenses will reduce 
the level of detail transfered to the image, but do not change the 
coherence properties of the partial waves. 
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3. Applications 



(corresponding to the L3 edge): 



3.1. Single atom 

For single, individual atoms, all final states are inde- 
pendent of one another and hence uncorrected. In addi- 
tion, in the absence of a (strong) external magnetic field, 
states with the same L, but different M or S can be con- 
sidered degenerate. For the XDOS, this means 



kn 
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In the case of no spin-polarization, the also do not 
depend on S and the sum over S, S' can be carried out 
and the spin-dependence of a can be dropped. Using the 
orthogonality relations for the Wigner 3j symbols, a short 
calculation then yields 



~ aa , =47r(2L + l)(2j + l) 



(18) 

Hence, for single, isolated atoms, S is diagonal, i.e., no 
correlations occur, even in the untransformed basis. Also 
note that all \i have the same weighting, i.e., any basis 
transform that only mixes the different /i of any given A 
but does not mix the A or L blocks leaves S unchanged. 
Owing to the incoherent summation over different \a) in 
eq. [l6j in the real space image this ultimately leads to 
terms of the forrrj^^ |e^| 2 = 1 [12^ which gives rise to 
circular intensity profiles regardless of the symmetries of 
the target's initial or final states. 

3.2. Energy-loss magnetic chiral dichroism 

Since its discovery in 2006 [6 , interest in the energy 
loss magnetic chiral dichroism (EMCD) technique has been 
growing steadily. Using EMCD, one can determine the 
magnetic properties of the sample [13] , similar to the X-ray 
magnetic circular dichroism which is a standard method in 
the synchrotron. The factorization approach outlined here 
can also be applied to EMCD. 

For the sake of simplicity, we assume here a fully spin- 
polarized (5 s i5 s ,i) dip ole- allowed (A = A' = 1) transition 
from an initial p (Z = 1) to a final d (L = V — 2) state, as is 
the dominant contribution to the L-edge in most common 
magnetic materials. In addition, we assume that states 
with same L, but different M are (mostly) degenerate, as 
in the isolated atom case. Hence, the XDOS reads 

^2 D LMS ( D I?M'S f ) = D 2 S L 2S L '2SmM'Ss^S s ,i. (19) 
kn 

Under these assumptions, S becomes a 3 x 3 matrix for 
both j = 1/2 (corresponding to the L2 edge) and j = 3/2 
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As in the single-atom case, S is already diagonal in the 
spherical harmonics basis. Here, however, different \i have 
different weights. This means that transforming to any 
other basis will introduce off-diagonal elements (only the 
identity matrix is invariant under rotations). Hence, the 
spherical harmonics basis is the only "physical" basis for 
EMCD. 

Moreover, the S matrices given above can be inter- 
preted as a homogeneous average signal on which the \i- 
dependent EMCD signal is superimposed: 
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This immediately shows two features common to EELS 
and EMCD. On the one hand, the homogeneous average 
signal exhibits the typical, statistical 1:2 intensity ratio of 
the L2:L3 edges. On the other hand, the EMCD signal 
(the absolute magnitude of which is independent of j in 
this simple case) shows the typical sign reversal between 
L2 and L3 edges. 

3.3. Crystals 

In crystals, the situation is more complicated and sim- 
ple toy-models are insufficient to grasp them completely. 
Hence, one needs sophisticated calculations to derive the 
XDOS that take into account the full crystal structure. 
Analogous to the examples given above (in particular to 
the isolated atom case), one can argue that for high-symmetry 
crystals (such as cubic systems), the XDOS may become 
diagonal in M, M' [HI [14] and hence many of the cross- 
terms vanish0 

Hence, we will use the Oxygen K-edge of Rutile (Ti02), 
a tetragonal system, as test case in this work. Fig. [TJl 
shows a schematic of the unit cell, while fig. shows the 
partial density of states (pDOS) for Oxygen as calculated 
by WIEN2k. From it, the lifting of the degeneracy of the 
different p orbit als is already evident. 

For this system, WIEN2k produces 89 non-negligibl^] 
XDOS components at Ep + 4 eV, whereas at Ep + 7 eV it 



5 Assuming elastic scattering effects are negligible. 



6 The complete investigation of the effects of crystal symmetries 
on the XDOS and S is beyond the scope of this work. 

7 Here, elements are considered non-negligible if they are larger 
than l%o of the largest element 
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Figure 1: (a) schematic of the unit cell for Rutile (gray: Ti, red: 
O). The lattice constants are a = 4.594 A and c = 2.958 A. The 
arrows show the symmetry-adapted local coordinate systems (red: x, 
green: y, blue: z). (b) pDOS of the Oxygen p-states as calculated by 
WIEN2k. The gray bars show energies used for simulations in this 
work. 



Figure 2: Real-space intensity of the exit wave after propagation 
of an incident plane wave through a one unit-cell thick crystal ori- 
ented in [0 1] zone axis at 200 kV acceleration voltage. Only dipole- 
allowed transitions were taken into account, (a) shows the image at 
an energy loss of Ep + 4 eV, whereas (b) shows the image at an en- 
ergy loss of Ep + 7 eV. The inset shows the projected unit cell with 
Ti atoms in gray and O atoms in blue. 
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Figure 3: Real-space intensity of the exit wave after propagation of 
an incident plane wave through a 10 nm thick crystal oriented in 
[0 1] zone axis at 200 kV acceleration voltage. A 24 mrad objective 
aperture was used. Only dipole-allowed transitions were taken into 
account. The image is taken at an energy loss of Ep + 4 eV. The 
inset shows the projected unit cell with Ti atoms in gray and O 
atoms in blue. 

produces 100 no n- negligible elements in the spin-unpolarized 
case. In the simplest case of taking into account only 
dipole-allowed transitions (A = X = 1), the 3x3 ma- 
trix S has 5 non- vanishing entries, which are reduced to 3 
after diagonalization. 

Fig. [2^i shows the simulated exit wave function inten- 
sities (corresponding to an ideal lens system) for a single 
unit cell after an energy loss of E F + 4 eV. The Oxygen p y 
orbitals are clearly visible to be pointing in the directions 
of the green axes in fig. |TJl. 

Likewise, fig. [2)3 shows the simulated exit wave function 
intensities at an energy loss of Ep + 7 eV. There, the 
Oxygen p z orbitals are clearly visible. They naturally are 
rotated by 90° with respect to the p y orbitals and are 
pointing in the direction of the blue axes in fig. [TJi. 

Fig. [3] shows the situation for a 10 nm thick crystal 
and an objective aperture of 24 mrad. It demonstrates 
that these results are not only of theoretical interest, but 
should be measurable in real instruments. The contrast 
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can be estimated to be 96 % (after subtraction of the pre- 
edge background, e.g., using the three- window method). 

More importantly, non-dipole transitions can easily be 
taken into account as well. For the Oxygen K-edge, the 
most relevant non-dipole transition is the monopole tran- 
sition from the Is state to final states with s symmetr^] 
Here, the number of non- negligible terms was reduced from 
10 to 4. Fig. [4] compares the intensities and images of 
monopole- allowed transitions, dipole-allowed transitions, 
and the coupling term between the two. A similar effect 
was predicted recently for X-ray absorption spectrometry 

Interestingly, the coupling term gives intensity varia- 
tions of almost ±20 % of the dipole-allowed transitions 
and is much larger than the monopole transitions them- 
selves. This acts in a way very similar to s-p hybridization, 
yielding an asymmetric image. 

4. Conclusion and Outlook 

In this work, we demonstrated a method to factorize 
both the MDFF and the rMDFF by means of a matrix di- 
agonalization. This was shown to yield obvious numerical 
advantages by reducing the number of terms to include in 
image calculations. Moreover, the diagonalization leads to 
a new set of basis vectors that are helpful to elucidate the 
physics underlying the scattering process. 

The new factorization method was applied to the iso- 
lated atom case, EMCD, and a Rutile crystal to show 
its versatility. In particular, the isolated atom and the 
EMCD cases could be treated analytically, giving results 
from which important properties such as the L2iL3 ratio 
or the sign reversal of the EMCD effect could be seen im- 
mediately. 

For the Rutile crystal, it was shown that with latest- 
generation TEMs, it should be possible to directly map 
atomic orbitals, e.g., using energy filtered TEM (EFTEM) 
with high spatial resolution. However, contrary to the 
common assumption that non-dipole transitions are unim- 
portant, it was shown that particularly the monopole- 
dipole coupling term can change the measured signal sig- 
nificantly. 

This new technique gives rise to exciting new possi- 
bilities. For example, it could be used to directly study 
the electronic structure of defects, interfaces, or other low- 
symmetry objects. 
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Figure 4: (a) - (d) Real-space intensities of the exit wave after propa- 
gation of an incident plane wave through a one unit cell thick crystal 
oriented in [0 1] zone axis at 200 kV acceleration voltage, (a) shows 
only monopole contributions (contrast-enhanced by a factor of 25), 
(b) shows only dipole contributions (contrast-enhanced by a factor 
of 1.2), (c) shows the total intensity, and (d) shows the coupling con- 
tribution (contrast enhanced by a factor of 7; red indicates positive 
values while cyan represents negative ones). The images are taken 
at an energy loss of Ep + 5.2 eV. The inset shows the projected unit 
cell with Ti atoms in gray and O atoms in blue, (e) Traces of the 
different contributions. The places of the trace are marked by yellow 
lines 



8 The pDOS for d states as produced by WIEN2k that would be 
accessible by quadrupole-allowed transitions is negligibly small. 
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Appendix A. Hermiticity of S 

For S to be hermitian, the equation 

must hold. 
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